General Update

The Sea Level Rise pressure layer was updated with new data for the year 2013. Previous data went through 2012. In addition, the rescaling of this pressure was recalculated to use the 99.99th percentile of raw data as the maximum. All values at or above the 99.99th percentile are given a value of 1.

Procedure

  1. Create raster of original data (came as .nc)
  2. project to mollweide
  3. Calculate cumulative sea level rise over all years
  4. Clip all negative values (values that indicate decreasing sea level)
  5. log transform
  6. resample to 1km
  7. rescale using 99.99 percentile
  8. Interpolate and replace NA cells with interpolated values (use Python - arcpy for this)
  9. Clip out ocean using ocean raster at 1km cell size
  10. Create raster of just interpolated cells

Data

Data for mean sea level rise (in mm) between January 1993 and June 2014 was downloaded on 1/12/15 from AVISO.

1. Create raster from original NetCDF file

# The following 3 lines were done once to create raster from downloaded NetCDF. 
# No longer need to run these

    #library(ncdf4)

    #r <- raster('raw/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.nc')
    #writeRaster(r,'tmp/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.tif')


# Read in raw data

    r <- raster('tmp/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.tif')

    plot(r,col=cols,main='Mean Annual Sea Level Rise (mm)\n1993-2014')

2. Reproject raster to mollweide and rotate to extent -180 to 180 degrees

  r = rotate(r) #rotate raster: -180 to 180
  
  #define initial projection. Helps avoid errors when reprojecting to mollweide
  projection(r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" 

  mollCRS <- CRS('+proj=moll')

  slr_moll <- projectRaster(r, crs=mollCRS, over=T)  

  plot(slr_moll,main='Sea Level Rise - Mollweide',col=cols)

3. Multiply mean annual sea level rise by number of years in dataset to get a total change in mm between 1993-2014

# (3) Multiply annual rate by all years in dataset

    # data used to create r comes from January 1993 - June 2014 (21.5 years)
    # Multiply values by 21+5/12 = 21.41667

      slr_moll <- slr_moll*21.41667

      plot(slr_moll,col=cols,main='Aggregate sea level rise (mm)\n1993-2014')

      histogram(slr_moll,main ='Distribution of values')

4. All values below zero indicate a decrease in sea level rise.

We choose to disregard these as the impact of sea level rise is focused only on the increase in sea level

# (4) Clip all negative values to 0

      slr_moll[slr_moll<0]<-0

      plot(slr_moll,main='SLR - negative values set equal to 0',col=cols)

      histogram(slr_moll,main ='Distribution of data after removing negative values')

5. Log transform the layer

    slr_moll_log = calc(slr_moll,fun=function(x){log(x+1)})

    plot(slr_moll_log,main='Sea level rise \nLog Transformed',col=cols)

    histogram(slr_moll_log,main='Distribution of data after log transform')

6. Resample layer to 1km resolution

Original resolution of the data is 0.25 degrees. To resample the data to 1km resolution, a template ocean raster that is at the desired resolution (~1km) and contains all cells included within OHI regions, is used.

    #ocean is a raster with all land clipped out - at 1km with value of 1
    ocean = raster(file.path(dir_N,'model/GL-NCEAS-Halpern2008/tmp/ocean.tif'))
    plot(ocean,main='Ocean Raster at 1km resolution \nUsed for resampling and clipping')

    slr_1km = resample(slr_moll_log,ocean,method='ngb')

    plot(slr_1km,main='SLR resampled to 1km' ,col=cols)

7. Rescale using the reference point (99.99 quantile)

All values at or above the 99.99th quantile (6.24183381230408) are set to equal 1.

    #get reference point
   
    ref = quantile(slr_1km,prob=0.9999)

    #normalize by the reference point - cap all values greater than 1 to 1
    r_resc <- calc(slr_1km,fun=function(x){ifelse(x>ref,1,x/ref)})

    plot(r_resc,main='Rescaled Sea Level Rise Pressure',col=cols)

8. Interpolation

Interpolating the data was done using ArcGIS. All NA cells were interpolated using nearest neighbor. The python script can be found here

  r_int = raster('tmp/slr_moll_log_1km_rescaled_int.tif')
  plot(r_int,main='Sea level rise interpolation with nearest neighbor',col=cols)

9. Clip out ocean & Finalize layer

Using the same ocean template raster as in step 6, all land is clipped out after interpolation.

  r_final = mask(r_int,ocean)

  plot(r_final,main='Sea Level Rise \n Final layer',col=cols)

10. Create interpolated cells raster

The cells that were filled in using the nearest neighbor interpolation are identified in a single raster.

interp = mask(r_final,slr_1km,inverse=TRUE)
plot(interp,col=cols,main='1km cells interpolated using Nearest Neighbor')